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. We performed Monte Carlo simulations to calculate the Berezinskii-Kosterlitz-Thouless (BKT) 

' temperature Tbkt for the two-dimensional planar rotator model in the presence of nonmagnetic 

impurity concentration (p). As expected, our calculation shows that the BKT temperature decreases 
as the spin vacancies increase. There is a critical dilution pc ~ 0.3 at which Tbkt ~ 0. The effective 
interaction between a vortex-antivortex pair and a static nonmagnetic impurity is studied analyt- 
ically. A simple phenomenological argument based on the pair-impurity interaction is proposed to 
justify the simulations. 
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PACS numbers: 75.30.Hx, 75.40.Mg, 75.10.Hk, 74.76.-w 



I. INTRODUCTION 



> 
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I The planar rotator (PR) model in two dimensions is a prototype for several physical systems as for example high 
■ temperature superconductors and granular superconductors. The PR model supports topological excitations and 
' although there is no long range order at any finite temperature, it undergoes a BKT phase transition driven by 
"Y^i the unbinding of vortex-antivortex pairs. In short the BKT picture of the phase transition is as follows. At low 
temperature spin waves are the relevant excitations of the system. Spin-spin correlation functions fall off slowly with 
distance, free vortices do not exist but pairs strongly binded. Vortices pairs can not disorder the system significantly 
since they affect only close spins. As the temperature is rised, the distance between vortex-ativortex pairs grows until 
' O I Tbkt- Then free vortices exist, the system is disordered and the spin-spin correlation function falls off exponentially. 
^ ' The hamiltonian describing the model is 

O 



- - ^ Ji^jSi.Sj, (1) 



' where i and j enumerate sites in a square lattice, Jij is an exchange coupling and Si = [Sf, Sf} = \S\ {cos0i,sm9i} 
5^ , is a two dimensional spin vector. Of course, Hamiltonian{|l|) describes an ideal system, in which each site of a regular 
square lattice is occupied by a spin vector S. However, impurity and/or defects are present in any material sample. 
In fact, the effect of impurities on superconductors has been of theoretical and experimental interest in its own right 
for a long time. Particularly, the interaction of topological excitations with spatial inhomogeneities is of considerable 
importance from both theoretical and applied points of view. Fot,gxample, solitons near a nonmagnetic impurity in 
2D antiferromagnets cause observable effects in EPR experimentflu. In this scenario it would be important to study 
the effects of the presence of nonmagnetic sites diluted in magnetic materials. In a recent work, M61, Pereira and 
PiresH, have studied the interaction between a static spin vacancy and a planar vortex and they have shown that the 
effective potential experienced between the two defects is repulsive. It indicates that the presence of spinless atoms 
on the magnetic plane may affect the BKT critical temperature. The main goal in this paper is to consider the effect 
of magnetic dilution to the BKT temperature by using numerical and analytical methods. To take into account the 
presence of nonmagnetic impurities in our model (Eq |^) we can replace some spin vector S" by a = creating a 
vacancy at that lattice site. First we consider that the spin vacancies are randomly distributed on the sites of the 
lattice. The case in which the spin vacancies are grouped into a cluster will also be analyzed in order to compare 
with the random case. 



The paper is organized as follows: in section II, we describe the model and the Monte Carlo (MC) method. In 
section /// we present the MC results. In section IV, the continuum theory is used to study the vortex-pair-impurity 
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interaction and a simple heuristic argument to justify the MC results is presented and section V contains a summary 
and final comments. 



II. BACKGROUND 



We consider in this work a quenched site diluted PR model. In order to introduce dilution we define a variable 
ai with the following properties: It is 1 if site i is magnetic and otherwise. To acomodate this changes we have to 
modify Eq |l| as 



H = —J ^ ^ aiSi.UjSj — — J ^ ^ (JiXTj 



(2) 



<ij> 



<ij> 



The precise determination of the BKT temperature is a difficult task due to absence o£|Sharp peaks in the thermo- 
dynamic quantities. One way to extract Tbkt was suggested by Weber and MinnhagenQcl by calculating the helicity 
modulus defined as 



T = 

where F is the free energy and A is a small twist across the system in one direction. Using Eq || we get 

1 „ 1 



T 



2{N - n) ^ hT{N - n) 



smi 



1,3 



(3) 



(4) 



where N is the volume of the system, n is the number of non-magnetic sites, e^j is the vector pointing frompsite j to 
site i and i; is a unit vector pointing along the x-direction. The Kosterlitz renormalization-group equationgj lead to 
the prediction that T jumps from the value (^) Tc to zero at the critical temperature, 



lim 

T^Te kbT 



(5) 



To calculate the quantity T we use a Monte Carlo (MC) approach using a standard Metropolis algorithm with periodic 
boundary conditionsQ. In order to reach the thermodynamical equilibrium we performed long runs of size 100 x Lx L, 
where L is the linear size of the lattice. The temperature was varied in steps of size AT = O.IK. Each point in our 
simulations is the result of the average over 2 x 10^ independent configurations. In the figures showing the results of 
our simulations, when not indicated, the error bars are smaller than the symbols. 

Figure 1 shows the results from MC simulations of T for lattices with 5% of impurities and sizes L = 30, 60 and 
80. The straight line represents (;|) T. The crossing point between this line and T gives an estimate of the BKT 
temperature. Of course, this estimate becomes more accurate as the lattice size increases. However, as we can see in 
figure 1, the lattice of size L = 60 gives already a good result adequated for our purposes. From now on we use the 
following, fcf, = 1 and the symbol Tbkt is used for Tc{p = 0), i.e, Tbkt = Tc{p = 0). 



III. MONTE CARLO RESULTS 



In this section we present the results obtained by MC simulations. First, we distribute the nonmagnetic impurities at 
random in the lattice sites. Figure ^ contains the helicity modulus as a function of the temperature considering several 
values of the impurity concentration {p). It is also shown the straight line representing the function (^) T. As noticed 
before the intersection of this line with the value of each T, gives T^ for the corresponding impurity concentration. We 
observe that Tic{p) decreases with increasing p. Since the helicity modulus is a measurement of the phase correlations 
of the systemO, it is not surprising that these correlations are strongly affected by the dilution. It can be understood 
as follow: if we remove a spin from the lattice, the nearest neighbors of that spin will have coordination number 
of three, one less than in the bulk. The spins in the boundary have larger fiuctuations than the spins in the bulk 
lowering the spin correlations. We should expect that the fluctuation becomes appreciable disordering the system for 
large enough nonmagnetic concentrations up to a critical value where the BKT temperature goes down to zero. In 
figure 3 we show the BKT temperature as a function of the nonmagnetic impurity concentration. Note the abrupt 
fall of the critical temperature for pc ~ 0-3. We also performed MC simulations for the case in which the nonmagnetic 
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FIG. 1: Helicity modulus T as a function of temperature for lattices with sizes 30x30, 60x60,80x80 and with 5% of nonmagnetic 
impurities randomly distributed. The solid line is the curve (-) T and the dashed lines are only guides to the eyes. 
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FIG. 2: Helicity modulus T as a function of temperature for lattices size 60x60 with 0%, 5%, 10%, 15%, 20%, 25%, 27%, 28%, 29% 
and 30% of nonmagnetic impurities randomly distributed. The solid line is the line (§) T and dashed curves are guides to the 
eye. 



impurities are clustered for p = 0.2 and 0.3 (see figure ||). Note that in this case, the critical temperature practically 
does not depend on the impurity concentration (Tc(0.2) = Tc(0.3)). In fact this is an expected result. Since the 
nonmagnetic cluster is confined in a region of size p x and the boundary grows as p x L, meaning that spins are 
still strongly correlated driving the BKT transition even for large values of p. A comparison between the two cases 
is shown in figure ^Note the considerable difference between them. Due to the short range of the spin interactions, 
only the spins near the boundary of the cluster will become infiuenced by the vacancies and hence the correlations of 
the rest of the system will have a behavior almost independent of the vacancies. It must not affect considerably the 
vortices that are formed far way from the cluster and the phase transition occurs normally. 



IV. VORTEX-ANTIVORTEX-IMPURITY INTERACTION 



In this section we discuss the effect of nonmagnetic sites on the vortex-antivortex structure. The interaction between 
the topological excitation and a single nonmagnetic impurity below the critical temperature may help us to understand 
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FIG. 3: The BKT transition temperature behavior as a function of nonmagnetic impurity concentration, based on the MC 
simulations results showed in figure 2. The dashed curves are guides to the eye. 
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FIG. 4: Helicity modulus T as a function of temperature, for lattices with 20% and 30% of nonmagnetic impurities grouped 
in a cluster, compared with the helicity modulus results for lattice with 20% of nonmagnetic impurities randomly distributed 
and lattice without impurities. The solid line is the line (-^^ T and dashed curves are guides to the eye. 



in more detail the phase transition mechanism. In the continuum limit, Hamiltonian (|^) can be written as 

He = (1/2)J / {yefd?x. (6) 



Following referenced, to take into account the absence of one spin in the lattice site we modify as 

Hi = (1/2)J / {VefV{r)cPx , (7) 



where V{r) is a localized potential given by: V{f) — \ ii \ f — fo\ > a, and V{r) = if |r — rd < a. Here, the 
nonmagnetic site is placed at Tq and a stands for the lattice constant. This lack of magnetic interaction inside the 
circle of radius a, means that a spin located at fo was removed from the lattice. The equation of motion obtained 
from (|^) is 

V{r)'^'^e = -yV{r).Vd. (8) 

In polar coordinates, the vectors r and are written as (r, (j)) and (r,,, (/)o) respectively. Then, the gradient of the 
potential is 

VF(r) = a[fcos{a — \4> — 0o|) + <f)sin{a — |0 — 4>o\)]5{r — fo — a), (9) 



FIG. 5: A vortex pair configuration with "center of mass" iocated at tlie origin and size R = 6a and a spin vacancy located at 
(1,3). The pair configuration is deformed for large distance if vortex centers are near the impurity, increasing considerably the 
system energy. 



where S is the Dirac delta function and a is the angle that the vector a, with origin at the point fo and end at a point 
on the circumference of the potential (|a| = a), makes with the vector Tq. In the limit a ^ 0, we write 



\7V{f) « a[fcos{a) + (j)sin{a)]d{f — Tq 



(10) 



where cos(a) and sin(a) are anisotropic coupling constants. A vortex-antivortex pair solution with "center of mass" 
at the origin is given by = arctan[{y — P)/x] — arctan[{y + P)/x], where R = 2P is the distance between the 
vortex centers. The energy of a pair is i?2i/ = tt^ J + 2TTjln{R/a). Note that the energy i?2i/ increases with increasing 
R implying an attractive force between vortices of opposite sign. Suppposing 9i is the deformation introduced in 
the vortex-antivortex structure by the absence of a spin at fo we can write 9 = 62,^ + Oi- Rewriting 621, in polar 
coordinates and substituting 9 = 92v + 9i into Eq. (0) , we obtain 



In 



''0 



where 



Pa[rrQC0s{2(j)Q — 0) + rP^coscj) — TqCos^o ~ -P^?'ocos(^o] 
|f - fo|(p4 + ^4 + 2P^4cos{2(j>Q)) 



(11) 



(12) 



The configuration of this deformed vortex-antivortex pair is shown in figure ^. As the vortex (or antivortex) center 
approaches the nonmagnetic impurity, the pair structure becomes more and more deformed, indicating that there is 
a repulsive interaction potential between each vortex core and the spin vacancy. This phaenomenum is in agreement 
with the results of Reffl, where the calculations took in consideration a single vortex. In order to understand how 
is the effective interaction potential between the two defects, we substitute 9 — 6*21/ -I- 9i into Eq.(0) to calculate the 
energy of the pair- impurity system Epi. Unfortunately, the integral in Eq.(0) can not be done analytically for a 
general impurity position, but in the special case the spin vacancy is located at the center of mass we can solve it 
exactly. Using the dominant terms, the effective potential is given by 









And' 












a 



(13) 



where Vef / = Epi — E2u , and d is the lattice size. Since P is the distance between the vortexj-(or antivortex) center and 
the spin vacancy, this expression is very alike with the effective potential obtained in RefJ3, between a single vortex 
and a nonmagnetic impurity. Note that the effective interaction potential increases with decreasing R — 2P, implying 
a repulsive force between vortices and impurities. In fact, the spin vacancy force obtained from Eq. ([T3|) acts as a 
"repulsive force" weakening the coupling strength between the bound vortices, and becomes stronger as P decreases. 
Here, the nonmagnetic impurity must repel simultaneously the two vortices in a pair, affecting the spin field for large 
distances (see fig. 5). For a lattice of size d, the effective potential ( p^ is a minimum only if P — > d/2{R d), 
showing the tendency of a complete separation of the vortices in a pair due to the presence of the vacancy. We 
conclude that static spin vacancies repel vortices, independently if they are free or bound into pairs. Based on the 
above results, we propose a phenomenological model to explain the behavior of the BKT temperature as a function 
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of the impurity concentration. As discussed above a nonmagnetic site can induce a a repulsive potential between a 
pair vortex-antivortex in such a way we can have the two scenario. If the nonmagnetic impurity is in between the pair 
vortex-antivortex the effective repulsive potential created tends to unbind the pair. On the other hand, if the impurity 
is not in between the pair the force in the nearest vortex will be stronger than in the other and the tendency is to 
increase the vortex-antivortex attraction leading to the annihilation of the pair. Then, impurities may induce cither 
vortex-antivortex unbinding process or pair annihilation. In a system containing a random distribution of impurities 
one can expect a lower density of vortices at any temperature than in a pure system due to the annihilation of pairs. 
Beside that, the unbinding of vortices-antivortices should occurs at lower temperature inducing the BKT transition. 

Hence, we may expect a critical nonmagnetic impurity concentration in which vortex pairs are not more formed and 
the BKT critical temperature goes to zero. The situation is different for the case in which the nonmagnetic impurities 
are clustered. In this case, vortex pairs will be excited far way from the cluster in order to minimize their energies and 
the cluster would have only a small influence on the vortex-antivortex unbinding. The critical temperature should 
not be much affected. The results presented in figure H confirms this conjecture. 



V. SUMMARY 



We have performed Monte Carlo simulation for the diluted planar rotator model in a square lattice. We have found 
that the BKT temperature decreases with increasing impurity concentration and that there is a critical impurity 
concentration pc = 0.3 at which the transition temperature goes to zero. The interaction between a vortex-pair and a 
static spin vacancy was studied in the continuum approximation. By considering the decoupling of vortex pairs induced 
by impurities we argumented that the BKT critical temperature should decrease, justifying the MC simulations. Our 
results may also have applications for granular superconducting films such as the ceramic high- Tc materials. These 
systems could be modeled as two-dimensional (2D) Josephson-junction arrays because such films contain large number 
of Josephson boundaries between the small superconducting grains forming a complex Josephson-junction network. 
However, the actual situation is not so ideal as the perfect array since grains with different sizes and orientations are 
arranged almost randomly. This makes the model with vacancies more realistic than the usual perfect array. The 
results can also be extrapolated to models with three spin components such as easy-plane and XY magnets. In these 
cases, the problem with impurities could be still more interesting since they have a true dynamics. In fact, it can shed 
some light over the important question about the origin of ^hf ..'"'^r'^i^"]. Flf "i^ iii the dynamical spin-spin correlation 
function in the two dimensional anisotropic Heisenberg modelElLE30't3aEj. However, much work has to be done in 
order to understand those effects. 
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